/*
 * GenClst_Control.cpp
 *
 *  Created on: May 28, 2011
 *      Author: nzhao
 */

#include "share/GenClst_Control.h"

int  * INIT_STATES;
int  * STATE_J_LIST;
int * NUM_CLST;
int  TOTALCLST;
int  ** CLST_INDEX;
Cluster_List CLUSTERS;
Sublist     *sblst, *splst;

void	generate_linkmap();
int		isLink(int i, int j);
void	clst_info_convert();
void	export_clst_info();

void generate_clusters()
{
	int i;
	if( CALCULATED_ORDER > 1)
		generate_linkmap();

	CLUSTERS.initialize();

	for(i = 1; i < CALCULATED_ORDER; i ++)
		CLUSTERS.generate_from_order(i);

	CLUSTERS.find_sub_clusters();

    clst_info_convert();
    CLUSTERS.clear();

//    export_clst_info();
}

void generate_linkmap()
{
	int i, j, count, a;
	LINK_MAP = new int [NBATH * ( NBATH + 1 ) / 2 ];

	count = 0; a = 0;
	for ( i = 0; i < NBATH; i++)
	{
		for (j = 0 ; j <= i; j++)
		{
			LINK_MAP[count] = isLink(i, j);
			a +=LINK_MAP[count];
			count++;
		}
	}
	cout << "Mean linke number = " << 2.0 * a / NBATH << endl;
}

int isLink(int i, int j)
{
	double dist2, dx, dy, dz;
	int res;

	if(i==j)
	{
		return 0;
	}
	else
	{
		dx =  COORDINATES[ 3 * i + 0] - COORDINATES[ 3 * j + 0];
		dy =  COORDINATES[ 3 * i + 1] - COORDINATES[ 3 * j + 1];
		dz =  COORDINATES[ 3 * i + 2] - COORDINATES[ 3 * j + 2];
		dist2  =  dx*dx + dy*dy + dz*dz;
		res = dist2 < THRESHOLDVALUE * THRESHOLDVALUE;

		if(IsMediated)
			res = res || ( MAX_MEDIATED_INT[i * NBATH + j] > MEDIATED_INT_THRESHOLDVALUE );

		return (res);
	}
}

void clst_info_convert()
{
    int i, j, q, clst_count = 0, len_sub, len_sup;

    NUM_CLST	= new int	[CALCULATED_ORDER];
    CLST_INDEX	= new int * [CALCULATED_ORDER];

    TOTALCLST = CLUSTERS.Total_Clst();
    for (i = 0; i < CALCULATED_ORDER ; i ++)
    {
        NUM_CLST[i ] = CLUSTERS.NUM_order(i);
        CLST_INDEX[i] = new int [ NUM_CLST[i] * ( i + 1 ) ];

        for (j = 0; j < NUM_CLST[i]; j++)
        {
                for (q = 0; q < i + 1; q++ )
                {
                        CLST_INDEX[i][j*(i+1) + q] = (CLUSTERS.index(clst_count)->mySpinIdx())[q];
                }
                clst_count ++;
        }
    }
    cout << "Cluster index is read" << endl;

    sblst = new Sublist[TOTALCLST];
    splst = new Sublist[TOTALCLST];

    for(i = 0; i < TOTALCLST ; i++)
    {
        len_sub = CLUSTERS.index(i)->my_subNum();
        len_sup = CLUSTERS.index(i)->my_supNum();

        sblst[i].sublen = len_sub;
        splst[i].sublen = len_sup;

        sblst[i].sublist = new int [len_sub];
        splst[i].sublist = new int [len_sup];
        for (j = 0; j < len_sub; j++)
        	sblst[i].sublist[j] = CLUSTERS.index(i)->my_Subs()[j];
        for (j = 0; j < len_sup; j++)
        	splst[i].sublist[j] = CLUSTERS.index(i)->my_Sups()[j];

    }
    cout << "SUB_IDX / SUP_IDX is read." << endl;

}

void export_clst_info()
{
	int i, len_sub, len_sup;
	char Cluster_file[200];
	strcpy(Cluster_file, WORKSHOP);
	strcat(Cluster_file, DEBUG_PATH);
	strcat(Cluster_file, "/clusters/clstinfo.dat");

	ofstream clstInfo(Cluster_file, ios::binary);
	if(!clstInfo) assert(0);

	clstInfo.write((const char *)&TOTALCLST, sizeof(int));
	clstInfo.write((const char *)&CALCULATED_ORDER, sizeof(int));
	clstInfo.write((const char *)NUM_CLST , CALCULATED_ORDER * sizeof(int));
	for(i = 0; i <CALCULATED_ORDER; i++)
	{
		clstInfo.write((const char *)CLST_INDEX[i] , (NUM_CLST[i] * ( i + 1 )) * sizeof(int));
	}
    for(i = 0; i < TOTALCLST ; i++)
    {
    	len_sub = sblst[i].sublen;
        len_sup = splst[i].sublen;

        clstInfo.write((const char *)&len_sub, sizeof(int));
        clstInfo.write((const char *)&len_sup, sizeof(int));

        clstInfo.write((const char *)( sblst[i].sublist ), len_sub * sizeof(int));
        clstInfo.write((const char *)( splst[i].sublist ), len_sup * sizeof(int));
    }


	clstInfo.close();

}

void generate_initial_state(int seed)
{
    int i, j, q;

//////////////////////////////////////////////////////////////////////////
// initialize a random number generator
    const gsl_rng_type * T;
    gsl_rng * r;

    gsl_rng_env_setup();

    T = gsl_rng_default;
    r = gsl_rng_alloc (T);

    gsl_rng_set(r, seed);

    INIT_STATES = new int [NBATH];

    double u;
    int pol = 0;
    for(i = 0; i < NBATH; i++)
    {
        u = gsl_rng_uniform(r);
        if( u < 0.5 )
            INIT_STATES[i] = 1; //spin up
        else
            INIT_STATES[i] = -1; //spin down

        pol += INIT_STATES[i];
    }
    cout << "polarization = " << pol << " / " << NBATH << " = " << (100.0 * pol) / NBATH << "%"<< endl;;


    STATE_J_LIST = new int [TOTALCLST];

    int base, st;

    int clst_count = 0;
    for (i = 0; i < CALCULATED_ORDER ; i ++)
    {
        for (j = 0; j < NUM_CLST[i]; j++)
        {
            base = 1;
            STATE_J_LIST[clst_count] = 0;
            for (q = i; q >= 0; q-- )
            {
                st = (1 - INIT_STATES [ CLST_INDEX[i][j*(i+1) + q] ] ) / 2;    //st: INIT_STATES[i]==1 --->0; -1---->1 by st = (1-init)/2

                    STATE_J_LIST[clst_count] +=  st * base;
                    base *= 2;
            }
            clst_count ++;
        }
    }

    gsl_rng_free (r);
}
